
rm(list=ls())

library(tidyr)
library(tibble)
library(stargazer)
library(dplyr)
library(haven)
library(ggplot2)
library(zoo)

#########################################################################
# REPLICATION FILES: MILITARIZATION AND PERCEPTIONS OF LAW ENFORCEMENT  #
# JOURNAL: BJPS                                                         #
# FIGURES IN SUPP. APP. FOR                                             #
# 1. DIFFERENCE IN AMCEs BY SUBGORUP                                    #
#########################################################################

#### STEP 1: LOAD RESULTS                     ####

# AGE

age_amcedif <-  read.csv("AMCE_Age.csv")

age_amcedif$upper95 <-age_amcedif$estimate + 1.96*age_amcedif$std.error
age_amcedif$lower95 <-age_amcedif$estimate - 1.96*age_amcedif$std.error

age_amcedif$upper90 <-age_amcedif$estimate + 1.645*age_amcedif$std.error
age_amcedif$lower90 <-age_amcedif$estimate - 1.645*age_amcedif$std.error  

# EDUCATION

edu_amcedif <-  read.csv("AMCE_Education.csv")

edu_amcedif$upper95 <-edu_amcedif$estimate + 1.96*edu_amcedif$std.error
edu_amcedif$lower95 <-edu_amcedif$estimate - 1.96*edu_amcedif$std.error

edu_amcedif$upper90 <-edu_amcedif$estimate + 1.645*edu_amcedif$std.error
edu_amcedif$lower90 <-edu_amcedif$estimate - 1.645*edu_amcedif$std.error  

# SEX   

gen_amcedif <-  read.csv("AMCE_Gender.csv")

gen_amcedif$upper95 <-gen_amcedif$estimate + 1.96*gen_amcedif$std.error
gen_amcedif$lower95 <-gen_amcedif$estimate - 1.96*gen_amcedif$std.error

gen_amcedif$upper90 <-gen_amcedif$estimate + 1.645*gen_amcedif$std.error
gen_amcedif$lower90 <-gen_amcedif$estimate - 1.645*gen_amcedif$std.error  


# HOMICIDES 5Y

hom_amcedif <-  read.csv("AMCE_Homicides.csv")

hom_amcedif$upper95 <-hom_amcedif$estimate + 1.96*hom_amcedif$std.error
hom_amcedif$lower95 <-hom_amcedif$estimate - 1.96*hom_amcedif$std.error

hom_amcedif$upper90 <-hom_amcedif$estimate + 1.645*hom_amcedif$std.error
hom_amcedif$lower90 <-hom_amcedif$estimate - 1.645*hom_amcedif$std.error

# IDEOLOGY   

id_amcedif <-  read.csv("AMCE_Ideology.csv")

id_amcedif$upper95 <-id_amcedif$estimate + 1.96*id_amcedif$std.error
id_amcedif$lower95 <-id_amcedif$estimate - 1.96*id_amcedif$std.error

id_amcedif$upper90 <-id_amcedif$estimate + 1.645*id_amcedif$std.error
id_amcedif$lower90 <-id_amcedif$estimate - 1.645*id_amcedif$std.error   

# INCOME

inc_amcedif <-  read.csv("AMCE_Income.csv")

inc_amcedif$upper95 <-inc_amcedif$estimate + 1.96*inc_amcedif$std.error
inc_amcedif$lower95 <-inc_amcedif$estimate - 1.96*inc_amcedif$std.error

inc_amcedif$upper90 <-inc_amcedif$estimate + 1.645*inc_amcedif$std.error
inc_amcedif$lower90 <-inc_amcedif$estimate - 1.645*inc_amcedif$std.error

# PERCEPTION

per_amcedif <-  read.csv("AMCE_Insecurity.csv")

per_amcedif$upper95 <-per_amcedif$estimate + 1.96*per_amcedif$std.error
per_amcedif$lower95 <-per_amcedif$estimate - 1.96*per_amcedif$std.error

per_amcedif$upper90 <-per_amcedif$estimate + 1.645*per_amcedif$std.error
per_amcedif$lower90 <-per_amcedif$estimate - 1.645*per_amcedif$std.error  

# TRUST

tru_amcedif <-  read.csv("AMCE_Trust.csv")

tru_amcedif$upper95 <-tru_amcedif$estimate + 1.96*tru_amcedif$std.error
tru_amcedif$lower95 <-tru_amcedif$estimate - 1.96*tru_amcedif$std.error

tru_amcedif$upper90 <-tru_amcedif$estimate + 1.645*tru_amcedif$std.error
tru_amcedif$lower90 <-tru_amcedif$estimate - 1.645*tru_amcedif$std.error    

# VICTIMS

vic_amcedif <-  read.csv("AMCE_Victim.csv")

vic_amcedif$upper95 <-vic_amcedif$estimate + 1.96*vic_amcedif$std.error
vic_amcedif$lower95 <-vic_amcedif$estimate - 1.96*vic_amcedif$std.error

vic_amcedif$upper90 <-vic_amcedif$estimate + 1.645*vic_amcedif$std.error
vic_amcedif$lower90 <-vic_amcedif$estimate - 1.645*vic_amcedif$std.error  

# PARTY ID

pid_amcedif <-  read.csv("AMCE_PID.csv")

pid_amcedif$upper95 <-pid_amcedif$estimate + 1.96*pid_amcedif$std.error
pid_amcedif$lower95 <-pid_amcedif$estimate - 1.96*pid_amcedif$std.error

pid_amcedif$upper90 <-pid_amcedif$estimate + 1.645*pid_amcedif$std.error
pid_amcedif$lower90 <-pid_amcedif$estimate - 1.645*pid_amcedif$std.error    


#### STEP 2: SET THEME AND LABELS                   ####

theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 8, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      axis.title.x =      element_text(size =11),
      legend.position = "none",
      panel.background = element_rect(fill = "#F5F5F5"),
      strip.background = element_rect(colour="black", fill="white")
    )
}

yylab1  <- c("Estimated Difference in AMCEs")


#### STEP 3: CREATE FIGURES                  ####


#---------------  FIGURE A.3 AGE ---------------------- #

age_amcedif$ylabel<-factor(age_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
age_amcedif$model<-factor(age_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))


age_amcedifplot= ggplot(age_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.1, 0.1)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.1, 0,0.1)) +
  facet_grid(.~ model) + 
  scale_x_discrete(name="")

plot(age_amcedifplot)

#pdf(paste("Age_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#age_amcedifplot = age_amcedifplot  + theme_bw1()
#print(age_amcedifplot)
#dev.off()

#ggsave("Age_AMCE_Differences.png", plot=age_amcedifplot, width=10,height=6, units = "in", dpi = 300)

#---------------  FIGURE A.4 SEX ---------------------- #

gen_amcedif$ylabel<-factor(gen_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
gen_amcedif$model<-factor(gen_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))


gen_amcedifplot= ggplot(gen_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.1, 0.1)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.1, 0,0.1)) +
  facet_grid(.~ model) + 
  scale_x_discrete(name="")

plot(gen_amcedifplot)

#pdf(paste("Figures/Choices for paper/AMCE Differences/Gen_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#gen_amcedifplot = gen_amcedifplot  + theme_bw1()
#print(gen_amcedifplot)
#dev.off()

#ggsave("Figures/Choices for paper/AMCE Differences/Gen_AMCE_Differences.png", plot=gen_amcedif, width=10,height=6, units = "in", dpi = 300)


#---------------  FIGURE A.5 EDUCATION ---------------------- #

edu_amcedif$ylabel<-factor(edu_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
edu_amcedif$model<-factor(edu_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))

edu_amcedifplot= ggplot(edu_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.1, 0.1)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.1, 0,0.1)) +
  facet_grid(.~ model) + 
  scale_x_discrete(name="")

plot(edu_amcedifplot)

#pdf(paste("Edu_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#edu_amcedifplot = edu_amcedifplot  + theme_bw1()
#print(edu_amcedifplot)
#dev.off()

#ggsave("Edu_AMCE_Differences.png", plot=edu_amcedifplot, width=10,height=6, units = "in", dpi = 300)

#---------------  FIGURE A.6 INCOME ---------------------- #

inc_amcedif$ylabel<-factor(inc_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
inc_amcedif$model<-factor(inc_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))


inc_amcedifplot= ggplot(inc_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.1, 0.1)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.1, 0,0.1)) +
  facet_grid(.~ model) + 
  scale_x_discrete(name="")

plot(inc_amcedifplot)

#pdf(paste("Inc_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#inc_amcedifplot = inc_amcedifplot  + theme_bw1()
#print(inc_amcedifplot)
#dev.off()

#ggsave("Inc_AMCE_Differences.png", plot=inc_amcedifplot, width=10,height=6, units = "in", dpi = 300)


#---------------  FIGURE A.7 TRUST ---------------------- #

tru_amcedif$ylabel<-factor(tru_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
tru_amcedif$model<-factor(tru_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))


tru_amcedifplot= ggplot(tru_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.1, 0.1)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.1, 0,0.1)) +
  facet_grid(.~ model) + 
  scale_x_discrete(name="")

plot(tru_amcedifplot)

#pdf(paste("Trust_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#tru_amcedifplot = tru_amcedifplot  + theme_bw1()
#print(tru_amcedifplot)
#dev.off()

#ggsave("Trust_AMCE_Differences.png", plot=tru_amcedifplot, width=10,height=6, units = "in", dpi = 300)

#---------------  FIGURE A.8 IDEOLOGY ---------------------- #

id_amcedif$ylabel<-factor(id_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
id_amcedif$model<-factor(id_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))


id_amcedifplot= ggplot(id_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.1, 0.1)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.1, 0,0.1)) +
  facet_grid(.~ model) + 
  scale_x_discrete(name="")

plot(id_amcedifplot)

#pdf(paste("Ideo_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#id_amcedifplot = id_amcedifplot  + theme_bw1()
#print(id_amcedifplot)
#dev.off()

#ggsave("Ideo_AMCE_Differences.png", plot=id_amcedifplot, width=10,height=6, units = "in", dpi = 300)

#---------------  FIGURE A.9 VICTIMIZATION ---------------------- #

vic_amcedif$ylabel<-factor(vic_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
vic_amcedif$model<-factor(vic_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))


vic_amcedifplot= ggplot(vic_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.1, 0.1)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.1, 0,0.1)) +
  facet_grid(.~ model) + 
  scale_x_discrete(name="")

plot(vic_amcedifplot)

#pdf(paste("Vic_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#vic_amcedifplot = vic_amcedifplot  + theme_bw1()
#print(vic_amcedifplot)
#dev.off()

#ggsave("Vic_AMCE_Differences.png", plot=vic_amcedifplot, width=10,height=6, units = "in", dpi = 300)


#---------------  FIGURE A.11 PERCEPTION ---------------------- #

per_amcedif$ylabel<-factor(per_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
per_amcedif$model<-factor(per_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))


per_amcedifplot= ggplot(per_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.1, 0.1)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.1, 0,0.1)) +
  facet_grid(.~ model) + 
  scale_x_discrete(name="")

plot(per_amcedifplot)

#pdf(paste("Perc_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#per_amcedifplot = per_amcedifplot  + theme_bw1()
#print(per_amcedifplot)
#dev.off()

#ggsave("Perc_AMCE_Differences.png", plot=per_amcedifplot, width=10,height=6, units = "in", dpi = 300)

#---------------  FIGURE A.12 HOMICIDES ---------------------- #

hom_amcedif$ylabel<-factor(hom_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
hom_amcedif$model<-factor(hom_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))


hom_amcedifplot= ggplot(hom_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.1, 0.1)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95),position="dodge",size=1, color="black") +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#969696", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.1, 0,0.1)) +
  facet_grid(.~ model) + 
  scale_x_discrete(name="")

plot(hom_amcedifplot)

#pdf(paste("Hom_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#hom_amcedifplot = hom_amcedifplot  + theme_bw1()
#print(hom_amcedifplot)
#dev.off()

#ggsave("Hom_AMCE_Differences.png", plot=hom_amcedifplot, width=10,height=6, units = "in", dpi = 300)

#---------------  FIGURE A.10 PID ---------------------- #

theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 8, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      axis.title.x =      element_text(size =11),
      legend.position = "bottom",
      panel.background = element_rect(fill = "#F5F5F5"),
      strip.background = element_rect(colour="black", fill="white")
    )
} 

pid_amcedif$ylabel<-factor(pid_amcedif$ylabel,levels=c("Light", "Skin Color:", "Male", "Gender:", "Assault Rifle", "Weapon:", "Military", "Uniform:"))
pid_amcedif$model<-factor(pid_amcedif$model,levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))

pid_amcedif$ylabel <- ifelse(pid_amcedif$party=="PRI",as.numeric(pid_amcedif$ylabel)-.3,pid_amcedif$ylabel)

pid_amcedifplot= ggplot(pid_amcedif, aes(y=estimate,x=ylabel)) +
  coord_flip(ylim = c(-.2, 0.2)) +
  geom_hline(yintercept = 0,size=.5,colour="darkgrey",linetype="solid") + 
  geom_linerange(aes(ymin=lower95,ymax=upper95, color=party),position="dodge",size=1) +
  geom_linerange(aes(ymin=lower90,ymax=upper90),size=1, color="#C0C0C0", linetype=1, position = position_nudge(x = -0.09)) +
  geom_point(size=3) +
  scale_y_continuous(name=yylab1, breaks=c( -0.2, 0,0.2)) +
  facet_grid(.~ model) + 
  scale_color_manual(name="Party ID:", values=c("#808080", "black")) +
  scale_x_continuous(name="",
                     breaks = 8:1,
                     labels = c("Uniform:",  "Military", "Weapon:", "Assault Rifle ",  "Gender:", "Male",  "Skin Color:", "Light"))


plot(pid_amcedifplot)

#pdf(paste("PID_AMCE_Differences.pdf",sep=""),width=7,height=5) 
#pid_amcedifplot = pid_amcedifplot  + theme_bw1()
#print(pid_amcedifplot)
#dev.off()

#ggsave("PID_AMCE_Differences.png", plot=pid_amcedifplot, width=10,height=6, units = "in", dpi = 300)
